using EcologicalNetworksDynamics
using Random, Plots, Distributions, DataFrames, StatsPlots
using EcologicalNetworksPlots, EcologicalNetworksThe previous tutorial focused on experiments where we manipulated the number of networks and various network parameters. This is one set of things we can change/vary in an in silico experiment. The other set of things we can change are features of the model, such as the shape of the functional response (see Tutorial 7), features of the environment such as the carrying capacity, or even empirical relationships that drive trophic structure and interaction strengths, such as the predator-prey mass ratio.
In this tutorial, we are going to implement three experiments. The first two will be ‘simple’ in that they vary only two things. The final example will implement a large experiment changing five features of the model.
You may want to start a new script in the project. We’ll need the following packages (they are already installed… so we just need using).
Experiment 1: Carrying Capacity and the Predator Prey Mass Ratio
Now we are set for our first experiment. Lets first establish the parameters we need to make the food web and do the experiment. We fix S at 20 and C at 0.15. We then create vectors of Z and K.
Z is the predator - prey mass ratio, and defines how much bigger or smaller the predators are from their prey. The data suggest it is between predators are between 10 and 100 times bigger than their prey see Brose et al 2006. This value interacts with setting trophic levels in the model.
The default setting for the models is 1 - i.e. all species are within the same order of magnitude, predators are not bigger than their prey. Here, we creat a vector of values to explore, from predators being smaller, to them being 10 or 100 x larger as the data suggests.
MORE ON HOW THE PPMR interacts with Trophic level?
#Fixed Parameters
S = 20
C = 0.15
# Variable Parameters
Z_levels = [0.1, 1, 10, 100]
K_levels = [0.1, 1, 10, 100]
# run this to get same results as in the document
Random.seed!(123)TaskLocalRNG()
Now, lets set up the collecting data frame.
df_collect = DataFrame(Z = [], K = [], FinalRichness = [], FinalBiomass = [], FinalStability = [])| Row | Z | K | FinalRichness | FinalBiomass | FinalStability |
|---|---|---|---|---|---|
| Any | Any | Any | Any | Any |
Now, set up the loop to use these variables and generate outputs
for z in Z_levels
for k in K_levels
println(" ***> This is iteration with Z = $z and K = $k\n")
fw = FoodWeb(nichemodel, S; C, Z = z)
Environment(fw, K = k)
B0 = rand(S)
params = ModelParameters(fw)
out = simulate(params, B0)
# collect info
fin_rich = foodweb_richness(out, last = 30)
fin_bio = total_biomass(out, last = 30)
stab = population_stability(out, last = 30)
push!(df_collect, [z, k, fin_rich, fin_bio, stab])
end
end ***> This is iteration with Z = 0.1 and K = 0.1
┌ Info: Species [16] went extinct at time t = 10.333220779245396.
└ 1 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 58.07510892697318.
└ 2 out of 20 species are extinct.
┌ Info: Species [8] went extinct at time t = 96.59439608251868.
└ 3 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 99.81800659688153.
└ 4 out of 20 species are extinct.
┌ Info: Species [10] went extinct at time t = 173.50017299410476.
└ 5 out of 20 species are extinct.
***> This is iteration with Z = 0.1 and K = 1.0
┌ Info: Species [7, 18] went extinct at time t = 7.79864262878258.
└ 2 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 12.302588972757128.
└ 3 out of 20 species are extinct.
┌ Info: Species [14, 19] went extinct at time t = 54.65909534164244.
└ 5 out of 20 species are extinct.
┌ Info: Species [10] went extinct at time t = 60.53358149873246.
└ 6 out of 20 species are extinct.
***> This is iteration with Z = 0.1 and K = 10.0
┌ Info: Species [5] went extinct at time t = 252.5627381565363.
└ 7 out of 20 species are extinct.
┌ Info: Species [16] went extinct at time t = 255.6762152446527.
└ 8 out of 20 species are extinct.
┌ Info: Species [18] went extinct at time t = 5.452845413338652.
└ 1 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 7.6238324278018.
└ 2 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 67.45898739065261.
└ 3 out of 20 species are extinct.
┌ Info: Species [19] went extinct at time t = 70.55144729932225.
└ 4 out of 20 species are extinct.
┌ Info: Species [11] went extinct at time t = 76.65337831298842.
└ 5 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 79.93906826954334.
└ 6 out of 20 species are extinct.
┌ Info: Species [9] went extinct at time t = 206.60040098574206.
└ 7 out of 20 species are extinct.
***> This is iteration with Z = 0.1 and K = 100.0
┌ Info: Species [9] went extinct at time t = 10.740211650168495.
└ 1 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 21.53600773237581.
└ 2 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 26.526137967578318.
└ 3 out of 20 species are extinct.
┌ Info: Species [10, 17] went extinct at time t = 29.65218241510892.
└ 5 out of 20 species are extinct.
┌ Info: Species [11] went extinct at time t = 40.93973683032451.
└ 6 out of 20 species are extinct.
┌ Info: Species [7, 3] went extinct at time t = 44.403037167488925.
└ 8 out of 20 species are extinct.
┌ Info: Species [5, 12] went extinct at time t = 47.422100875122986.
└ 10 out of 20 species are extinct.
┌ Info: Species [15] went extinct at time t = 53.521946406912726.
└ 11 out of 20 species are extinct.
┌ Info: Species [13] went extinct at time t = 88.23965733559986.
└ 12 out of 20 species are extinct.
***> This is iteration with Z = 1.0 and K = 0.1
┌ Info: Species [7, 11] went extinct at time t = 33.40288048145466.
└ 2 out of 20 species are extinct.
┌ Info: Species [16] went extinct at time t = 43.58846848478901.
└ 3 out of 20 species are extinct.
***> This is iteration with Z = 1.0 and K = 1.0
┌ Info: Species [13] went extinct at time t = 71.26258343675524.
└ 4 out of 20 species are extinct.
┌ Info: Species [9] went extinct at time t = 76.2910085126206.
└ 5 out of 20 species are extinct.
┌ Info: Species [17] went extinct at time t = 81.47041929142708.
└ 6 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 250.51901692767748.
└ 7 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 30.025171952720154.
└ 1 out of 20 species are extinct.
┌ Info: Species [16, 20, 17, 19] went extinct at time t = 41.4267587744762.
└ 5 out of 20 species are extinct.
┌ Info: Species [4] went extinct at time t = 61.788208267960805.
└ 6 out of 20 species are extinct.
┌ Info: Species [11, 18] went extinct at time t = 108.48651859473986.
└ 8 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 146.3080674035236.
└ 9 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 296.7203892892406.
└ 10 out of 20 species are extinct.
***> This is iteration with Z = 1.0 and K = 10.0
┌ Info: Species [7] went extinct at time t = 31.97167042395919.
└ 1 out of 20 species are extinct.
┌ Info: Species [20, 19] went extinct at time t = 44.704826238230446.
└ 3 out of 20 species are extinct.
***> This is iteration with Z = 1.0 and K = 100.0
┌ Info: Species [16] went extinct at time t = 351.1728007870926.
└ 4 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 93.92059343902136.
└ 1 out of 20 species are extinct.
┌ Info: Species [13] went extinct at time t = 232.64065602391264.
└ 2 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 257.55666477477075.
└ 3 out of 20 species are extinct.
***> This is iteration with Z = 10.0 and K = 0.1
***> This is iteration with Z = 10.0 and K = 1.0
┌ Info: Species [19] went extinct at time t = 176.22073341629988.
└ 1 out of 20 species are extinct.
┌ Info: Species [14] went extinct at time t = 447.84231449641817.
└ 2 out of 20 species are extinct.
***> This is iteration with Z = 10.0 and K = 10.0
┌ Info: Species [18] went extinct at time t = 257.2149928441204.
└ 1 out of 20 species are extinct.
┌ Info: Species [12] went extinct at time t = 362.2028229313007.
└ 2 out of 20 species are extinct.
***> This is iteration with Z = 10.0 and K = 100.0
┌ Info: Species [12] went extinct at time t = 392.3491322903875.
└ 1 out of 20 species are extinct.
┌ Info: Species [3] went extinct at time t = 108.36651973920425.
└ 1 out of 20 species are extinct.
┌ Info: Species [5, 4] went extinct at time t = 142.19180110075308.
└ 3 out of 20 species are extinct.
┌ Info: Species [11, 9] went extinct at time t = 152.16900575546077.
└ 5 out of 20 species are extinct.
┌ Info: Species [6] went extinct at time t = 162.39218592381178.
└ 6 out of 20 species are extinct.
┌ Info: Species [15, 10] went extinct at time t = 184.60725874208512.
└ 8 out of 20 species are extinct.
┌ Info: Species [7] went extinct at time t = 196.60065563305216.
└ 9 out of 20 species are extinct.
***> This is iteration with Z = 100.0 and K = 0.1
┌ Info: Species [16, 13] went extinct at time t = 254.46108282792977.
└ 11 out of 20 species are extinct.
┌ Info: Species [20] went extinct at time t = 324.694736964769.
└ 12 out of 20 species are extinct.
┌ Info: Species [19] went extinct at time t = 360.7907924164594.
└ 13 out of 20 species are extinct.
***> This is iteration with Z = 100.0 and K = 1.0
***> This is iteration with Z = 100.0 and K = 10.0
***> This is iteration with Z = 100.0 and K = 100.0
┌ Info: Species [4] went extinct at time t = 303.7164417912157.
└ 1 out of 20 species are extinct.
Wonderful. Now we are in a position to learn about two new plotting methods. First, let’s look at the data frame we’ve created.
df_collect| Row | Z | K | FinalRichness | FinalBiomass | FinalStability |
|---|---|---|---|---|---|
| Any | Any | Any | Any | Any | |
| 1 | 0.1 | 0.1 | 15.0 | 2.69924 | -0.0076855 |
| 2 | 0.1 | 1.0 | 12.0 | 2.13677 | -0.0011011 |
| 3 | 0.1 | 10.0 | 13.0 | 2.18971 | -0.000546612 |
| 4 | 0.1 | 100.0 | 8.0 | 1.50359 | -0.00038262 |
| 5 | 1.0 | 0.1 | 14.4 | 2.22238 | -0.136409 |
| 6 | 1.0 | 1.0 | 12.1333 | 2.64512 | -0.178709 |
| 7 | 1.0 | 10.0 | 16.8 | 2.67383 | -0.259915 |
| 8 | 1.0 | 100.0 | 18.7 | 3.54345 | -0.0968901 |
| 9 | 10.0 | 0.1 | 19.4 | 4.81281 | -0.232271 |
| 10 | 10.0 | 1.0 | 19.5 | 8.08706 | -0.101531 |
| 11 | 10.0 | 10.0 | 19.2667 | 4.24443 | -0.0289865 |
| 12 | 10.0 | 100.0 | 12.6667 | 4.70276 | -0.0880609 |
| 13 | 100.0 | 0.1 | 20.0 | 10.8534 | -0.42463 |
| 14 | 100.0 | 1.0 | 20.0 | 41.1869 | -0.421736 |
| 15 | 100.0 | 10.0 | 19.7 | 14.6763 | -0.750017 |
| 16 | 100.0 | 100.0 | 20.0 | 18.7396 | -0.433297 |
Visualising the experiment
One option here is to plot one of our Final Objects as the response variable against the valuse of Z and K. In R, we’d use ggplot2. Here we’ll use StatsPlots as we learned about in Tutorial 5. Can you make this work in the regular Plots syntax?
Let’s first look at a single plot of stability
@df df_collect plot(:K, [:FinalStability], group = :Z,
ylabel = "Stabilty",
xlabel = "Karrying Kapacity",
seriestype = [:scatter, :line],
legend = false)p1 = @df df_collect plot(:K, [:FinalStability], group = :Z,
legend = :bottomright,
ylabel = "Stabilty",
xlabel = "Karrying Kapacity",
seriestype = [:scatter, :line])
p2 = @df df_collect plot(:K, [:FinalBiomass], group = :Z,
legend = :bottomright,
ylabel = "Stabilty",
xlabel = "Karrying Kapacity",
seriestype = [:scatter, :line])
p3 = @df df_collect plot(:K, [:FinalRichness], group = :Z,
legend = :bottomright,
ylabel = "Stabilty",
xlabel = "Karrying Kapacity",
seriestype = [:scatter, :line])
# create a layout of 3 graphs stacked on top of each other.
plot(p1, p2, p3, layout=(3,1), legend = false)